Solvent free model for self-assembling fluid bilayer membranes: 
Stabilization of the fluid phase based on broad attractive tail potentials 



in 
o 
o 

(N 

Oh' 
(D 
00 

OO 



o 



o 



> 

OO 

0^ 
O 

o 



O 

o 



Ira R. Cooke and Markus Deserno 
Max-Planck-InsUtut fiir Polymerforschung, Ackermannweg 10, 55128 Mainz, Germany 

(Dated: February 2, 2008) 

We present a simple and highly adaptable method for simulating coarse-grained lipid membranes 
without explicit solvent. Lipids are represented by one head-bead and two tail-beads, with the 
interaction between tails being of key importance in stabilizing the fluid phase. Two such tail-tail 
potentials were tested, with the important feature in both cases being a variable range of attraction. 
We examined phase diagrams of this range versus temperature for both functional forms of the 
tail-tail attraction and found that a certain threshold attractive width was required to stabilize the 
fluid phase. Within the fluid phase region we find that material properties such as area per lipid, 
orientational order, diffusion constant, inter-leaflet flip-flop rate and bilayer stiffness all depend 
strongly and monotonically on the attractive width. For three particular values of the potential 
width we investigate the transition between gel and fluid phases via heating or cooling and find that 
this transition is discontinuous with considerable hysteresis. We also investigated the stretching of a 
bilayer to eventually form a pore and found excellent agreement with a recently published analytic 
theory. 

PACS numbers: 61.20.Ja, 81.16.Dn, 82.70.Uv 



I. INTRODUCTION 

Lipid bilayers are among the most versatile of nature's 
biomaterials. As the interface between the cell and its 
environment, or between organelles and the cytosol, they 
provide regulated transport of substances as small as pro- 
tons to as large as entire cells. Such a wide range of func- 
tional length scales is naturally studied via an equally 
broad range of methods. For example, at the smallest 
scale detailed quantum atomistic simulations are neces- 
sary to study the transport of ions or water across the 
membrane interface [1, 2], whereas at the opposite end 
of the length scale spectrum, analytic theory [3] and dy- 
namically triangulated lattice simulations [4] have been 
used to determine the shape behavior of whole vesicles 
under various conditions of pressure, volume and area, or 
even under hydrodynamic flow [5] . Between these two ex- 
tremes are problems at the so called "meso-scale" which 
include membrane fusion and rupture [6-8], domain for- 
mation in multi-component membranes [9-12], the cou- 
pling of membrane composition with curvature, or the 
interaction of membranes with colloidal or viral particles 
[13]. Such problems occur at relatively large length and 
timescales but at the same time require a particle based 
approach that reproduces the basic bilayer structure of 
the membrane. The combination of these two require- 
ments necessitates the use of coarse grained simulation 
approaches, in which groups of atoms are represented by 
single particles. Such coarse grained approaches vary in 
their level of detail from a single bead per lipid [14] up 
to quite detailed lipids with on the order of ten beads 
[15]. Naturally such a range of levels of coarse grain- 
ing goes along with a tradeoff between computational 
efficiency and level of detail. In this respect, the single 
most important determinant of model speed appears to 
be the presence or absence of explicit solvent. Unlike the 



two-dimensional membrane, the solvent is the hulk phase 
which fills the entire simulation box, and integrating its 
degrees of freedom can easily amount to more than 90% 
of simulation time. Naturally, one must include solvent 
when its effects are of inherent interest to the physics of 
the problem. In many cases however, its task is merely to 
mediate the hydrophobic attraction between lipid tails, 
and as such it is secondary to the overall purpose of the 
simulation. 

If one could simulate bilayer membranes without the 
need for explicit solvent, a vast increase in accessible 
length and time scales would result, yet despite the long 
and successful history of solvent free models in polymer 
physics, this approach has not yet been widely adopted 
for lipid bilayer simulations. This is because the mem- 
brane case displays one additional complication: Unlike 
polymers, whose structure is at the outset determined 
by chemistry, lipids first have to physically self-assemble 
into a two-dimensional fluid bilayer. This aggregation 
results from a balance between lipid entropy and the en- 
ergy of cohesion. Since in the solvent-free case the lat- 
ter stems from effective attractions (for instance between 
lipid tails), a physically meaningful balance will pose re- 
strictions on the interaction potentials. Indeed, the col- 
lective experience from the past has shown that simple 
choices {e.g. Lennard- Jones, LJ) do not lead to a fluid 
bilayer phase but only to "solid" bilayers at low temper- 
ature and low density ("gas") phases at high tempera- 
ture. Concluding that simple pair potentials are insufh- 
cient, researchers have then turned to the use of density 
dependent (multibody) interactions [6, 13, 14, 16, 17], 
angular dependent potentials [18] or highly tuned sets of 
Lennard-Jones like potentials [19] to stabilize the fluid 
bilayer phase without solvent (see Brannigan et al. for 
a recent review [20]). Unfortunately, each of these ap- 
proaches suffers from one or more significant drawbacks. 
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For example, the multibody approach introduces compH- 
cations for interpretation and measurement of thermo- 
dynamic quantities, while neither the angular dependent 
approach nor the use of tuned LJ potentials has so far 
led to bilayers for which unassisted self-assembly has been 
demonstrated. In addition to these technical problems, 
it is also notable that the reported bending stiffnesses 
are generally outside the experimentally reported range 
[21-23], being restricted to either relatively low ([14, 17] 
< bknT) or high ([18, 19] > 50A:bT) values. Thus, there 
remains a clear need for an efficient solvent-free bilayer 
model that does not suffer from such technical drawbacks 
and is also highly tuneable. 

Recently, two new solvent free models have appeared 
which set out to solve these problems. One has been pro- 
posed by Brannigan et al. [24] , the other one by us [25] . 
Both models display a wide parameter range in which 
simple pair attractions drive an unassisted self-assembly 
into a fluid phase, within which the bending stiffness can 
be easily tuned. They differ from all previous solvent 
free models in their use of attractive potentials that ex- 
tend somewhat further than a simple LJ-potential. These 
act either between all tail beads [25] or are restricted to 
special interface beads between hydrophilic head and hy- 
drophobic tail [24]. As we shall show in Sec. Ill, our 
experience with broad attractions of various functional 
forms strongly suggests that it is this feature which ul- 
timately enables a fluid bilayer phase for these strongly 
coarse grained systems. 

Although similar in spirit, these two models differ in a 
number of details. For instance, Brannigan et al. opt for a 
more detailed representation of lipid molecules compared 
to us, with the intention to capture local features of the 
bilayer stress profile more accurately, but at a concomi- 
tant price in efficiency. It therefore depends on the phys- 
ical problem under study, as well as on ones position in 
the detail vs efficiency tradeoff, which model is preferable 
in any given situation. A judicious choice then requires 
good knowledge of the physical properties of these mod- 
els, which have so far been outlined only rather briefly, 
and the present paper goes toward flUing this gap for the 
model developed by us. In Sec. II and III we briefly de- 
scribe our model and discuss its self-assembly properties. 
In particular, we support our claim that the long-ranged 
nature of the attractions are the feature of key impor- 
tance by showing that the qualitative physical proper- 
ties are robust against change of the specific functional 
form of the attraction. In Sec. IV a detailed account of 
the properties of the fluid phase and their variation with 
attractive width and temperature is given along with a 
description of the gel-fluid transition. Finally, in Sec. V 
we study the stretching and rupture of a bilayer sheet 
and find near perfect agreement with a simple theoreti- 
cal model developed by Farago [19] and also Tolpekina 
et al. [26] as well as with experimental data. 



II. BASIC PRINCIPLES OF THE MODEL 

We describe a model surfactant system which is based 
on the simple idea that the solvent mediated interaction 
between lipid tails can be represented by an effective at- 
tractive potential of sufficiently broad range. Such a po- 
tential can meet the two demands of (i) providing enough 
cohesive energy to drive assembly of a two-dimensional 
aggregate, while (ii) permitting enough lateral freedom 
for the lipid constituents to remain in a fluid state. There 
are many ways in which one might construct a coarse 
grained model based around the central principle of such 
a broad attractive potential. Below, we shall present two 
alternatives which differ in the exact form of their tail 
attractions. In showing both of these we merely seek to 
demonstrate that it is not the precise functional form 
that matters but merely the attractive range. 

In this paper we present the most coarse grained ver- 
sion of our model, because this is the most useful in terms 
of length scales obtainable. Nonetheless, it is worth not- 
ing that the same principles can trivially be extended 
to include more detailed lipids with greater numbers of 
atoms should this be required. Of course, one should al- 
ways remember that by increasing the number of lipid 
beads, the ratio of simulation lengths to real lengths be- 
comes less favorable. The lipids we use arc represented 
by one "head" bead followed by two "tail" beads. Their 
size is fixed via a Weeks-Chandler-Andersen potential 

with rc = 2^/^&. We use e as our unit of energy. In order 
to ensure an effective cylindrical lipid shape we choose 
fchoad.hcad = fchcad,taii = 0.95 (T and 6taii,taii = CT , where cr 
is our unit of length. The three beads are linked by two 
FENE bonds 

V1,ond(r) = -ifcbond rl, log [1 - {r/r^f] , (2) 

with stiffness fcbond = 30 e/cr^ and divergence length 
Too = 1.5 (7. Lipids are straightened by a harmonic spring 
with rest length 4cr between head-bead and second tail- 
bead 

Vbend(r) = ^kbendir - ^Cff , (3) 

which corresponds in lowest order to a harmonic bending 
potential ^fctend^^^?^ for the angle tt — between the 
three beads. We fixed the bending stiffness at fcbend = 
lOe/cr^. 

The absence of explicit solvent molecules and the hy- 
drophobic effect they would give rise to is compensated 
by an attractive interaction between all tail beads. We 
compared two alternative potentials that account for this 
effect (see insets to Figs. 1 and 2 ). The first of these 

{-e , r < rc 

-e cos2 ^^^^^ , rc < r < rc + Wc (4) 
, r > rc -I- ti^c 
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describes an attractive potential with depth e which for 
r > Tc smoothly tapers to zero. In this case, tuning the 
decay range Wc proves to be the key to obtaining a fluid 
bilayer state. 

The second alternative is based on the familiar 
Lennard-Jones potential but extends its range simply by 
inserting a flat piece of length W[ at the minimum 

VflatLjW = (5) 

{— e , r < Tc + u>f 

, r > Wf+ Went 

where our key tuning parameter is now the width W{ 
of the flat region. The potential is cut-off beyond wi + 
Wcut, where Wcut = 2.5 cr is set to the usual value. Note 
that tuning W{ achieves a broad potential T^iat lj simply 
by introducing a flat region, whereas in the case of Vcos 
the tuning parameter Wc varies decay range and shape 
simultaneously. 

The above model is sufficiently simple to allow imple- 
mentation with a variety of molecular dynamics (MD) 
integration schemes or even Monte-Carlo. In this work 
we performed MD simulations with a Langevin thermo- 
stat to obtain the canonical ensemble [27] (time step 
St = 0.01 T and a friction constant T — r^^ in Lennard- 
Joncs units [55]). Constant volume simulations were per- 
formed using a cuboid box with sides = Ly, Lz sub- 
ject to periodic boundary conditions. If needed, constant 
tension conditions were also implemented via a modified 
Andersen barostat [28] allowing box resizing in x and y 
dimensions only (with a box friction Fbox = 2 x 10~^t~^ 
and box mass within the range Q — 10~^ . . . 10"^). By 
employing such a barostat rather than a simple Monte 
Carlo box length move we aimed to preserve the correct 
fluctuation behavior of our system. 

All simulations were performed using the ESPResSo 
program [29] . Lipid membrane specific analysis and setup 
was done using the mbtools package which is included as 
part of the main ESPResSo distribution. 



III. BILAYER STABILITY AND SELF 
ASSEMBLY 
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FIG. 1: Phase diagram resulting from Vcos cohesion (Eqn. 4) 
in the plane of potential width Wc and temperature at zero lat- 
eral tension. Each symbol corresponds to one simulation and 
identifies difi'erent bilayer phases: x: gel; •: fiuid, +: unsta- 
ble. Lines are merely guides to the eye. The inset shows the 
pair-potential between tail lipids (solid line) and the purely 
repulsive head- head and head-tail interaction (dashed line). 
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FIG. 2; Phase diagram resulting from Vflat lj cohesion 
(Eqn. 5) in the plane of potential width Wf and tempera- 
ture at zero lateral tension. The meaning of all symbols is 
the same as for figure (Fig. 1) 



In this section we shall map out the conditions under 
which a tensionless fluid bilayer state is stable. To do this 
we identified the fluid phase in two different ways (see 
Figs. 1 and 2 for the following) . In the first method, a 
box-spanning bilayer was pre-assembled from 1000 lipids, 
and its equilibration under zero lateral tension was at- 
tempted. This resulted in one of three possible outcomes: 
a stable fluid bilayer, a gel phase with strongly increased 
lipid order and much lower diffusion constant, or com- 
plete breakup of the system to form a "gas" phase. The 
transition between gel and liquid occurred over a narrow 
temperature range and will be explored in more detail 
in Sec. IV C. The gas phase was always identified as the 



point at which the imposition of zero tension conditions 
resulted in a divergence in the box length. 

Our second method for identifying the fluid bilayer 
phase was designed to ensure that we were not artifi- 
cially stabilizing the bilayer due to pre-assembly. In this 
case we conducted constant volume simulations starting 
from a random "gas" configuration. Under all conditions 
which previously gave stable tensionless membranes, a 
bilayer patch quickly self-assembled, which, at the cor- 
rect box size could zip up with its open ends to span the 
box (see Fig. 3 for such a sequence). If the box was too 
big, the patch either remained free, or (sometimes) closed 
upon itself to form a vesicle. For rather large values of the 
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FIG. 3: Self-assembly sequence for the bilayer system with 
1000 lipids in a cubic box of side length 25a. Lipid cohesion 
was set to Wc/a — 1.4 (Vflat lj) and temperature to ksT — e. 
A random gas of lipids quickly forms small clusters which 
slowly coarsen and eventually "zip up" to form a box-spanning 
bilayer sheet. The numbers indicate the MD time. 



width parameter Wc or W{ the line between fluid bilayer 
and the isotropic "gas" state becomes less distinct. For 
example, we observed that for Wc ^ 1.6 tr and Wf > 0.4 a, 
self assembly to box spanning bilayers (at constant vol- 
ume) could occur well above the evaporation boundary 
and that immediately below this boundary we observed 
rather indistinct bilayers, with particularly high flip-flop 
rates and low orientational order. Although we will see 
in Sec. IV A 2 that these relatively disordered states still 
show a definite bilayer structure, we should not be sur- 
prised if strange behaviors are found in this small region, 
and we would therefore caution against the physical in- 
terpretation of such results without additional checks for 
model artifacts. 

Using the methods mentioned above, we determined 
the phase diagrams for both types of attractive tail po- 
tentials, Eqns. (4) and (5), shown in Figs. 1 and 2, re- 
spectively. The most important point to note is that in 
both cases we see a fluid bilayer region that broadens sig- 
nificantly as potential width is increased. In this sense, 
both phase diagrams are remarkably similar given the 
strong differences in the nature of the functional forms 
and tuning parameters used to obtain them. At a more 
detailed level, one can see differences in the shape of the 
transition lines between the two models. This is likely 
due to the fact that the two width parameters work en- 
tirely differently. In the case of the cosine attraction, 
Eqn. (4), the attractive gradient is actually varied along 
with Wc, whereas for Eqn. (5), varying W{ leaves the at- 



tractive shape unchanged and merely adds a region of 
zero force close to the particle. 

Based on the relative similarity of our results for two 
such radically different potentials, our expectation is 
therefore that other potentials with a broad attractive 
width should also share the important property of ex- 
hibiting a stable fluid bilayer phase. This finds further 
support in the observation that the recent model of Bran- 
nigan et al. [24] also uses a broad 1 /r^ attraction. Of 
course, a direct comparison is difficult because they re- 
strict this attraction to the interface bead between lipid 
head and tail. Notwithstanding the physical motivation 
of this choice as being in accord with knowledge of the 
lateral stress profile, it would be interesting to check 
whether fluidity is ultimately insensitive to this detail 
and rests on the long range alone. 

Faced with these possibilities, as well as others we 
might think of for the tail attraction, the question is now 
which to choose. Since our potential attempts to capture 
the effects of solvent exclusion, lipid-lipid interactions, 
and the fact that our 3 bead "lipid" is a highly coarse 
grained representation of the real thing, it is difficult to 
guess what its functional form should look like. Instead 
we could attempt to differentiate between models on the 
basis of their emergent physical properties; however, it 
turns out that both are highly tuneable over a similar 
range. Therefore our primary considerations are practi- 
cal ones. In this regard the cosine attraction, Eqn. (4), 
is preferable since it acts over a slightly shorter range 
than the broadened LJ potential and is therefore faster 
to compute. The remainder of this paper will therefore 
focus on it alone. 

IV. PROPERTIES OF THE FLUID PHASE 

We have characterized bilayers in both fluid and gel 
phases by several observables, including the lipid orienta- 
tional order parameter, cross-bilayer density profiles and 
bending modulus, as well as the dynamical quantities dif- 
fusion constant and flip-flop-rate. In Sec. IV A we shall 
define these quantities and explain in detail how they are 
measured. In Sec. IV B the results are presented as cross 
sections at constant temperature, while Sec. IV C exam- 
ines the fluid-gel transition via cross sections at constant 
potential width. 

A. Observables 

1. Orientational order parameter 

Lipids in the fluid Lq phase are on average oriented 
parallel to the bilayer normal. The amount of alignment 
can be quantified by an orientational order parameter 
defined by 

^=i(3(a, •n)2-l), (6) 
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where is the unit vector along the axis of the i hpid, 
n is the average bilayer normal and angular brackets in- 
dicate an average over all lipids. A completely isotropic 
system will have S = whereas a fully ordered crystalline 
bilayer will have S = 1. 

2. Cross-bilayer density profile 

In a well-defined bilayer the lipid distribution perpen- 
dicular to the bilayer plane is very regular. In particu- 
lar, each of the constituent beads should occupy a well 
defined vertical distance from the bilayer midplane. To 
investigate this aspect of fluid bilayer structure we have 
calculated the number density of beads p{z) as a func- 
tion of vertical distance z from the local bilayer midplane. 
We used systems with 4000 lipids at constant zero ten- 
sion with a lateral box size of Lx = Ly k 50 a. Although 
such a large system provides good statistics, it also in- 
troduces the problem of dealing with undulations. We 
solved this by first assigning lipids to a 16 x 16 grid in 
the xy-plane and measuring the height z with respect to 
the average local height for that grid cell. Failure to do 
so will overestimate the width of the distributions signif- 
icantly. 

Calculations of p{z) for weakly coarse grained [30] or 
fully atomistic lipid bilayers [31] typically show differ- 
ences between the shape and width of p{z) for each func- 
tional group in the lipid. Of course, a strongly coarse 
grained model such as ours cannot hope to reproduce 
such subtleties. Our main concern is rather to ensure 
that the bilayer is not merely a loose agglomeration of 
lipids but that these lipids are oriented approximately 
vertically and that they do not strongly interdigitatc 
(although such interdigitated phases do occur naturally 
under certain circumstances [32, 33] they are not the 
"norm" for a fiuid bilayer). From Fig. 4 one can ver- 
ify that a well defined bilayer structure is indeed present 
for our system. Each of the three lipid beads shows a 
sharp peak about its average z. While the width of these 
peaks broadens upon approaching the liquid-gas bound- 
ary, their location is relatively stable. Upper and lower 
head beads are separated by a distance of approximately 
4.5 cr, while the infiection points of the summed density 
are separated by about 5cr. This agrees with our ex- 
pectations for vertically oriented lipids having only mi- 
nor interdigitation. Further evidence for distinct bilayer 
leafiets can be seen in the summed bead density which 
shows a slight peak for each of the terminal tail beads 
and a minimum at the bilayer midplane. 

In order to quantitatively examine trends in peak 
broadening and overlap throughout the fluid phase region 
of our phase diagram we define an overlap parameter ^ 
as follows: 

where i and j are labels for the bead type (i. e., head. 
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FIG. 4; Profile of the density p as a function of vertical 
distance z from the bilayer midplane for a system of 4000 
lipids at constant zero tension and with simulation parame- 
ters fcaT = l.le and Wc = 1.6. Plotted lines are bead densi- 
ties for head beads (dotted line), first tail beads (dot-dashed 
line), terminal tail beads (dashed line), and the sum of all 
beads (solid line). 

taili, tail2) and Vl = J2i I dz(pi(2))^ is a normalization 
factor. From this definition we see that complete overlap 
(or bead equivalence) is indicated by \E' = 1, whereas 
a rigid crystalline structure with no overlap would give 
^' = 0. 



3. Bending modulus 

One of the key material properties of a macroscopic 
membrane is its bending stiffness, which measures the 
energetic cost per unit area of imposing a local curvature. 
More precisely, the classical continuum description [34, 
35] in the limit of almost ffat membranes states that the 
energy of a deformed piece of membrane is given by 

E J dxdy [K{Ahf + j:{\/hf] , (8) 

where k and S are bending modulus and lateral tension, 
respectively, and where h{x,y) describes the height of the 
membrane above some reference plane ( "Monge gauge" ) . 
If we expand h{x,y) in Fourier modes according to 

h{r)=Y,hge"i-^ with q=^{nx,ny) (9) 
q 

and insert the result into Eqn. (8), we see that the en- 
ergy reduces to a sum of uncoupled harmonic oscillators 
with the modes hq as the degrees of freedom. From the 
equipartition theorem we then get the power spectrum 
of modes as [3] 
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FIG. 5: Asymptotic q ^ scaling of the power spectrum (|ftq|) 
for the bilayer system with Wc/a — 1.4 and ksT/e — 1.0. 



It is standard practice to obtain the bending modu- 
lus from a fit of the measured fluctuation spectrum to 
Eqn. (10). However, some care has to be taken here. 
First, Eqn. (8) is a continuum description and thus re- 
quires us to focus on large l ength scales. But for wave- 
vectors q smaller than Qc = the dominant influence 
is the tension S, and (|ft.q|) ~ is insensitive to k, so 
larger wave- vectors than qc are needed. However, once 
1/q becomes comparable to the bilayer thickness, simple 
continuum theory breaks down and further effects {e.g. 
protrusion modes [36]) set in. Identifying the character- 
istic scaling of the bending regime over a sufficiently 
wide range thus requires qc - and therefore the lateral 
tension S - to be as small as possible. Unfortunately, it 
turns out to be extremely hard to eliminate any remain- 
ing tension by adjusting the simulation box size by hand, 
since the tension depends very sensitively on bilayer area 
(see Sec. V below). We avoided this difficulty by instead 
using a modified Andersen barostat [28] to simulate in 
an ensemble of constant zero tension. Furthermore, to 
get away from microscopic lengths we took systems four 
times as big as the ones we used for mapping the phase 
diagram (4000 lipids, L ~ 50 a). Note that reaching the 
continuum limit in MD simulations is not trivial, since 
the relaxation time of bending modes scales as q~'^ . 

After setting up the bilayer, we first waited until ten- 
sion, box length, and energy had equilibrated (which 
took typically 10^ r for fluid systems). Then on the or- 
der of 100 configurations separated by 500 r were used 
to measure the mode spectrum. The bilayer mid-plane 
was identified by tracking the tail-beads and interpolat- 
ing their vertical position onto a 16 x 16 grid. Possible 
stray lipids had to be excluded from this procedure. A 
Fast Fourier Transform then yields the power spectrum 
but this requires one additional correction due to 
amplitude under-sampling on a grid which we briefiy dis- 
cuss in Appendix A. Fig. 5 provides a typical example of 
such a mode spectrum from which we can clearly see the 



asymptotic scaling, but also the deviations at large q. 
In this case length scales exceeding L « 20 cr {i. e., about 
four times the bilayer thickness) are required to reach the 
asymptotic regime. Hence, a simulation of smaller sys- 
tems (1000 lipids, L 25a) would not suffice to obtain a 
fluctuation spectrum which could in any meaningful way 
be fitted to Eqn. (10). 



4. Diffusion constant 

Calculating the in-plane diffusion constant Db for 
lipids in a solvent-free bilayer system is not as trivial 
a task as it might at first seem. For a strictly two- 
dimensional diffusive system we would of course have 



Dy. 



4^(["ll''^ 



it + At)-r\\4t)f) , (11) 



where r\\ i{t) is the projection [56] of the position vector 
of the lipid at time t into the bilayer plane. At is the 
time difference over which diffusion is probed, and the an- 
gular brackets indicate an average over all N lipids as well 
as configurations separated by a time difference of At. 
The difficulty, though, is that the system is really three- 
dimensional and that a small but nonvanishing fraction 
a ^ 1 of all lipids (typically a ~ 1 — 2 %) resides in the 
gas phase surrounding the bilayer. Even though they are 
the minority phase, they of course diffuse much faster and 
might thus contribute significantly to the average mean 
squared horizontal distance traveled. Unfortunately it is 
not straightforward to eliminate such stray lipids from 
the average {■ ■ ■)i over all lipids in Eqn. (11), since this 
would require us to check on the positions of all lipids 
during all times between t and t + At — information that 
is not necessarily available. 

Fortunately there is a way to check Eqn. (11) which 
does not require us to follow lipids, namely, by looking 
at the entire distribution function P{s) of lipid displace- 
ments s = (Ar||)^. If we think of lateral lipid diffusion as 
being due to two independent simple diffusion processes 
- one with bilayer diffusion constant Dy, and one with an 
effective lateral diffusion constant Dg ^ Z?b through the 
gas phase - a simple calculation would suggest 



P{s,At) = (1-a)- 



4DbAt 



,-s/4_DgAt 

4D^At 



(12) 



Indeed, a histogram of squared traveled distances shows 
exactly this double exponential decay, with the short dis- 
tance behavior dictated by the real bilayer diffusion con- 
stant _Db (see Fig. 6). We tested in different cases (in- 
cluding ones with a rather large fraction of stray lipids) 
that Dh is independent of At and in fact given by the 
diffusion constant obtained from the much simpler anal- 
ysis via Eqn. (11). In contrast, Dg depends essentially 
inversely on At, which rather than suggesting that diffu- 
sion through the gas phase somehow slows down reminds 
us that free diffusion through three-space is ultimately 
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FIG. 6: Scaled distribution of squared displacements s of 
lipid molecules for the system with Wc/cr = 1.6, ksT/e — 1.1 
and At = 2000 r (crosses), At = 6000 r (open circles), and 
At = 12000 r (filled circles). The lines are fits to Eqn. (12). 
These always yield a « 0.02 and values of the two diffusion 
constants as given in the inset. While Db is constant, the 
data are compatible with Dg approaching D^, with an 1/At 
asymptotics (inset, dashed line). 



cut short via readsorption in the bilayer or its periodic 
image. 

Having thus established that fast moving stray hpids 
have no significant influence on the bilayer diffusion con- 
stant as measured via Eqn. (11), we subsequently used 
this simpler analysis to obtain Db- In order to make 
best possible use of our data we calculated Db by taking 
the average over all possible values of the starting time 
t for each interval length At (which varied from 2000 r 
up to the entire length of an equilibrium simulation run 
IOOOOOt)) and then took the weighted average over 
all possible values of At. 



5. Flip-Flop rate 

A lipid molecule will not stay forever in the particular 
monolayer in which it presently resides; rather, there is 
a particular probability per unit time, the flip-flop-rate 
r, that it changes the monolayer. Let N+{t) and -/V-(t) 
be the number of lipids in the upper or lower monolayer 
at time t, which were present in the upper layer at time 
t = 0. These numbers satisfy the Master equations 

N±{t + dt) ^ N±{t){l-rdt) + rdtN^(t) , (13) 

The total number of such lipids is conserved, N+{t) + 
N_{t) = N/2, so we obtain N±{t) = -r[2N±{t) - N/2]. 
The solution satisfying the initial conditions iV+(0) = 
N/2 and iV_(0) = is N±{t) = ^N{1 ± e'^''*). The 
same considerations hold for lipids which at t — were 
in the lower monolayer. Hence, the total fraction f{t) of 
lipids which at time t reside in the same monolayer as 



they did at time i = is given by 

/W = ^(l + c-^'-*) . (14) 

This fraction is easily measurable in simulations. A fit 
to Eqn. (14) then yields the flip-flop-rate r. Notice that 
the probability density for flipping at time t is given by 
re~^*, and thus the average time between flip-flops is 
it) = 1/r. 

B. Constant temperature cuts 

Fig. 7 summarizes the values of orientational order pa- 
rameter S (Eqn. 6), overlap order parameter 4" (Eqn. 7), 
area per lipid a, and bending modulus k for constant 
temperature scans along the isotherms k^T / 1 = 0.6, 0.8, 
1.0, and 1.1. In order to more directly compare between 
results at different temperature we applied the rescaling 
w'^ = l + Wc — w]? where w]^ represents the value of Wc at 
the liquid-gas transition line for each respective temper- 
ature (we chose 1 as a reference value rather than to 
permit the construction of log-scale plots) . This allows 
us to scale out the most obvious effects of temperature 
and brings the results for all of the isotherms much closer 
to a common trend. 

We first note that the trends for all observables are 
monotonic, showing a consistent change from bilayers 
close to the liquid-gas boundary to those close to the 
liquid-gel boundary. Starting with the structural param- 
eters, S and 5", we see that both indicate an increase 
in order with larger potential range Wc. This is not im- 
mediately obvious, because broadening of the attractive 
potentials could also give the lipids more lateral freedom. 
Yet, the increase in S indicates that lipids fluctuate less 
around their average vertical position, and this goes hand 
in hand with a concomitant decrease in the overlap \E' be- 
tween the vertical density profiles for individual beads. 
In fact, it turns out that the average bead positions are 
roughly constant, hence the main mechanism by which 
'I' can increase is via the broadening of p(z) peaks for 
individual beads (see Fig. 4). Close to the liquid gas 
boundary the bilayer order is rather low. Indeed, visual 
inspection of bilayers with > 0.1 confirms that these 
are indeed very "fuzzy" . Although a clear bilayer struc- 
ture can still be seen in such cases, we would caution 
against using these in the attempt to model real lipid 
systems. 

The increase in order becomes more understandable 
when looking at the average area per lipid. Extending 
the range of the cohesive potential leads to an overall 
lateral contraction of the bilayer, thus explaining the re- 
duced fiuctuations and thus the behavior of S and '5. It 
is remarkable that the area per lipid of all four isotherms 
agrees after rescaling. This indicates that the lipid den- 
sity depends purely on the distance from the liquid-gas 
phase boundary rather than the absolute value of tem- 
perature or Wc individually. 
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FIG. 7: Basic static properties of the fluid bilayer phase as a function of w'^ and ksT, where w'^ indicates a rescaled attractive 
potential (w'^ = l + Wc — w]?). w^^ is the value of Wc on the liquid-unstable (gas) transition line. Each plot shows four isotherms; 
ksT = 0.6e (filled squares), ksT = 0.8e (asterisks), ksT = l.Oe (open circles), ksT = l.le (filled circles). The values of w]? for 
each of these isotherms were 0.185, —0.025, —0.2 and —0.27 respectively. In all cases statistical errors were smaller than the 
size of plotted points. 



The area a per lipid can also be used to map the coarse- 
grained length scale a to experimental lengths. For real 
phospholipid membranes values around 0.75 nm^ for the 
area per lipid are typical [31, 37], while our simulations 
give values in the range 1.1 — 1.5 a^. Assuming that one 
coarse-grained lipid is equivalent to one real lipid this 
gives a mapping of roughly a ~ 0.7 — 0.8 nm. An al- 
ternative mapping can be obtained by comparing a typ- 
ical bilayer thickness of roughly 5 nm with the measured 
width of the overall lipid density of approximately 5 a 
(see Sec. IV A 2 and Fig. 4), which gives a ~ Inm. Such 
good agreement between these two mappings indicates 
that our very simple 3 bead lipids are actually remark- 
ably close to the aspect ratio of real lipids. 

The observables presented so far reflect mostly local bi- 
layer properties. In contrast, the bending stiffness is an 
observable which, even though it ultimately derives from 
local bilayer properties, yields physics which is accessi- 
ble by large scale continuum calculations. Much of the 
theoretical modeling of fluid membranes hinges on the 



remarkable fact that on sufficiently large length scales 
they can be described by idealized surfaces with a very 
simple energy density, for which the bending stiffness k is 
in almost all cases the only relevant modulus [3, 34, 35]. 
Reproducing experimentally meaningful and easily tune- 
able values for this modulus is therefore one of the key 
requirements for any coarse-grained membrane modeling 
that aims at bridging the gap between local and global 
scales. 

Using the procedure described in Sec. IV A 3 we cal- 
culated K for each of the four isotherms k-gT/e = 0.6, 
0.8, 1.0, 1.1 and for all values of Wc in the fluid phase. 
The two most important conclusions from these data are 
the following: First, the range of accessible values for 
the bending stiffness coincides exactly with the experi- 
mentally interesting range for usual phospholipids [22]. 
Second, the value of k can be easily tuned via one pa- 
rameter, the potential range Wc- Beyond that, it is quite 
remarkable that plotting n/k^T against the rescaled po- 
tential width w'^ collapses all data points onto one mas- 
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FIG. 8: Diffusion constant as a function of rescaled potential 
width = 1 + Wc — uic^. Symbols and shifts w^^ are the same 
as in Fig. 7. 
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FIG. 9: Flip-flop-rate r as a function of rescaled potential 
width vj'^ — l + Wc —w]?. Symbols and shifts wi^ are the same 
as in Fig. 7. 



to the unusual shape of the phase diagram (see Fig. 1) 
between Wc/cr — 0.8 and Wc/a — 1.1. 

A typical value for the diffusion constant of lipids in 
real phospholipid membranes is about 1 /im^/s [38]. Tak- 
ing an average value of about 0.01 cr^/r from our data, 
and using the length mapping cr w 1 nm (see above) we 
obtain a time scale mapping of t « 10 ns. Although we 
do not place great quantitative weight on such an ap- 
proximate calculation, it does serve to illustrate that the 
timescales accessible by us are extremely long, being of 
the order of milliseconds. This is long enough, for ex- 
ample, to allow vesicles to self assemble or fuse, and for 
macroscopic phase separation to take place. 

In the case of flip-flop rates we found a strong de- 
pendence on temperature which is in accordance with 
the fact that this is a thermally activated process. The 
most important point to note however is that the range 
of actual values for rf obtained by us are many orders 
of magnitude faster than those typically found for artifi- 
cial phospholipid bilayers in experiments [39, 40]. Such 
a large discrepancy is not as alarming as it may at first 
seem since the rate rf is exponentially dependent on the 
activation energy for flip Ef (Arrhenius Law). We can 
therefore account for a very large difference in rf by a rel- 
atively small discrepancy in Ei. Nevertheless, it is clear 
that the lipids in our model undergo flip-flop too easily 
compared with real lipids. If one were specifically inter- 
ested in this aspect of the dynamics, this would clearly 
represent a problem, however in many other cases it pro- 
vides an advantage because the system will approach 
equilibrium more rapidly. If an accurate flip-flop rate 
is important we anticipate that our model lipids could 
be made to flip much less readily simply by increasing 
the chain length slightly (eg 4 bead lipids) and imposing 
a much stronger head-tail repulsion. 



C. Constant Wc profiles: Gel-Fluid Transition 



ter curve - with the notable exception of k^T/e = 0.6, 
which also stands out for the diffusion constant (see be- 
low). As a guide to the eye, a single line is shown for 
the three warmest isotherms and a separate line for the 
coldest one. Just as for the area per lipid it seems thus 
that the bending modulus depends only on the distance 
from the liquid-gas boundary. This is plausible, since the 
bending modulus is inversely proportional to the lateral 
compressibility, which again depends on the lipid density. 

Both dynamical properties, the bilayer diffusion con- 
stant _Db (Fig. 8) and the flip-flop-rate rf (Fig. 9) show 
a clear exponential decay with increasing w'^. Since we 
know that for free diffusion D oc T, a simple rescal- 
ing by temperature brings the diffusion constants for all 
isotherms into rough agreement, with the notable ex- 
ception of fceT/e = 0.6. We have seen that this same 
isotherm also gives an anomalous trend in the bilayer 
stiffness, and one can speculate that this may be due 



At constant values of Wc one can vary the tempera- 
ture to observe both liquid-gel and liquid-gas transitions. 
Since we studied lipid phase behavior at vanishing lat- 
eral tension, the liquid-gas boundary is necessarily sharp 
in our case. However, it is worth mentioning that un- 
der constant volume conditions alternative low density 
phases such as spherical or cylindrical micelles can be 
observed. Using the present model we have in fact ob- 
served such micellar phases, and together with one more 
control parameter - overall lipid density - their phase be- 
havior could be studied as well. This, however, is not the 
purpose of our present work and will be presented else- 
where. Here we concentrate on the transition between 
liquid and gel phases at vanishing lateral tension. There 
are several different ordered phases which collectively can 
be referred to as gels [32, 41], but we shall not attempt to 
identify all of these. Indeed, our rather simple three bead 
model was neither designed to reproduce such subtleties 
of lipid ordering nor would we actually expect to observe 
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FIG. 10: Variation of the area per lipid a across the gel-fluid 
phase boundary. Each figure shows a cooling-heating hystere- 
sis for a particular value of Wc- From top to bottom the values 
of Wc/cr used were 1.0, 1.4, and 1.6. Arrows indicate the direc- 
tion of temperature change. The rate of temperature change 
was 2.5x 10^® e/ksT for the top plot and 5x 10"'' e/ksr for the 
bottom two plots. The three vertical lines in the uppermost 
plot indicate the temperatures where the order parameter of 
Fig. 11 has been measured. 



the full zoo of ordered bilayer phases. 

In Fig. 10 we show the variation of lipid area a with 
temperature for three values of Wc- In all cases we chose 
to employ continuous temperature scans in order to allow 
the barostat to smoothly follow the substantial changes 
in box size involved. The important question then arises: 
"is our rate of cooling slow enough?" . We checked this 
by comparing plots of a us r for runs with four different 
cooling rates lO"-*, 10~^ 2.5 x 10"^, and 10~^ (in units of 
t/k-QT), from which one can see a clear deviation for the 
fastest rate but qualitatively very similar results for all 
slower ones. The only difference is that slightly sharper 




FIG. 11: Probability density p{z) of the height difference z be- 
tween a lipid and its 6 immediate neighbours. Solid, dashed, 
and dotted curves correspond to the temperatures indicated 
by lines 1, 2, and 3 in Fig. 10, respectively. 



transition boundaries can be seen as the rate is lowered, 
even among the slowest three rates. As a compromise 
between speed and resolution we chose the intermediate 
rate 5x e/k^T for runs with Wc/cr = 1.4 and Wc/ct = 
1.6 which display only a single transition. For runs with 
Wc/cr = 1.0, where two transitions must be resolved, we 
chose the slower rate, 2.5 x 10~^ e/ksT. 

Before focusing on the results at each value of Wc it 
is interesting to note that all runs show quite a strong 
hysteresis across the transition boundary as is typical 
of first order transitions. However, there is also a long 
tail to the hysteresis which appears during cooling. This 
is most likely due to the fact that the kinetics during 
gelling is strongly determined by the slow healing of de- 
fects. Indeed, we often observed such defects, many of 
which could be seen to dissappear during further slow 
cooling. 

Turning now to the results at Wc = 1.0 we find that 
there are two clear and sudden transitions during cooling 
as well as the reverse for heating. During cooling, both 
transitions involve a contraction of the area a and an in- 
crease in the orientational order S (not shown); however, 
if we look at the diffusion constant we find that it 
decreases suddenly from about 2 x 10""^ cr^/r to approx- 
imately 4 X 10~^<7^/t in the first (higher temperature) 
transition but does not decrease further during the sec- 
ond (lower temperature) transition. In fact, the very 
small value of the diffusion constant is numerically hard 
to determine accurately, but the overall drop by about 
two orders of magnitude is probably robust and corre- 
sponds well to what is known for typical phospholipid 
bilayers [38]. Thus we have a first transition from fluid 
to gel and a second transition at a lower temperature 
to a different gel phase. In order to understand what is 
actually occurring in the gel-gel transition we require a 
suitable order parameter, and it turns out that local lipid 
packing is what one has to look at. We calculated a his- 
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togram of height differences Az between a hpid head and 
that of its six nearest neighbors. The results are shown in 
Fig. 11 for bilayers at temperatures k^T/e — 0.38, 0.42, 
and 0.46 corresponding to each of the three phases. Both 
the fluid phase and the high temperature gel phase show 
a distribution peaked about Az = 0; however, there is a 
striking change as we go to the phase at fceT/e — 0.38 
which shows three peaks at Az = ±0.5 cr and 0. This 
shows that the lipids are now packed in an off-centered 
manner rather like oranges in a crate. Puzzlingly, how- 
ever, an equal proportion of lipids are found at Az ~ 
which would not occur if a perfect packing was obtained. 
This might be due to the presence of a still large number 
of defects. 

At larger values of Wc there is just a single transition 
from the fluid to the gel phase. This gel phase is packed 
similarly to the low temperature gel phase for Wc — 1.0 
since it also exhibits a three peak structure in the his- 
togram of Az. 

Finally, the plots of area per lipid vs temperature 
also contain information besides the evident phase tran- 
sitions. Their slope yields the value of the lateral thermal 
area expansivity, defined by 



For instance, the system with Wc/cr = 1.6 has a slope 
d{a)/dT K, Q.lQQa^ /{e/k-o) in the fluid phase and an 
area per lipid of 1.142(7^ at k-QT/e — 1. Identifying this 
temperature with room temperature, we obtain a thermal 
expansivity of ar ~ 2.2 x K^^ , which coincides 
remarkably well with the value measured by Kwok and 
Evans for fluid lecithin bilayers (2.4 x lO-^if-^) [42]. 

V. MEMBRANE STRETCHING AND PORE 
OPENING 

So far we have studied membranes under zero lateral 
tension. If we now apply extensional stress, the area per 
lipid will increase, up to the point where structural sta- 
bility of the bilayer is breached. Beyond a critical stress 
a pore can be nucleated, which then grows indefinitely, 
i e., the bilayer ruptures. This scenario has been de- 
scribed theoretically by Litster [43]. Important bilayer 
properties (e. g. the line tension of an open edge) could 
be extracted from observables such as the critical stress 
or the critical pore size, but it is evidently experimen- 
tally very difficult to study membranes at the brink of 
rupture. However, in a simulation it is very easy to per- 
form measurements in a different ensemble, namely one 
of constant area of the entire bilayer. Beyond some crit- 
ical strain one would now expect a pore to open, but 
since this now relieves much of the stress, pore growth 
stops and one obtains a stable pore of well-defined size. 
Assuming a harmonic extensibility of the bilayer itself, 
as well as a constant line tension at a pore rim, Farago 
[19] Tolpekina et al. [26] gave a simple theoretical model 
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FIG. 12: Bilayer tension E as a function of (projected) area 
Atot for a flat membrane sheet with Wc/a = 1.6 at ksT/e = 
1.1. The bold solid line is a fit to the model of Tolpekina 
et al. [26] (see also Eqn. (B5) in Appendix B) ; the fine solid 
and dashed curves indicate metastable and unstable branches, 
respectively. 



which relates the resulting pore size as well as the stress- 
strain relation to key bilayer properties such as the ex- 
tensibility and the line tension. We summarize the key 
results in Appendix B. 

Using our simulation model with the parameters 
Wc/<J — 1.6 and k^T/t — 1.1 we determined the equi- 
librium lateral tension S as a function of total box area 
{A — Lx 'X- Ly). For each value of the area this was done 
by placing 4000 lipids in a bilayer configuration spanning 
the xy-plane of the simulation box. After allowing for 
an equilibration time of 3 x 10** we then simulated for 
a further 3 x 10^, during which time the lateral tension 
S = —Lz{pxx +Pyy)/'^ was measured. The resulting plot 
of A S is shown in Fig. 12. Three main regimes are 
clearly distinguishable, separated by two values of the 
bilayer area. First, at a particular area Aq the tension 
vanishes. Boxes with a smaller area yield a negative ten- 
sion, i. e., a positive pressure. The bilayer is under com- 
pressional strain, which it evades very soon by buckling. 
Conversely, boxes with an area larger than Aq subject 
the bilayer to extensional strain and create a proportional 
increase in tension. At some particular value for the bi- 
layer area the energy stored in the extension, which grows 
quadratically with the strain, must exceed that of a bi- 
layer with a pore, since the pore size cannot grow faster 
than the strain, and the line tension is proportional to 
the square root of pore size. Once the pore opens it re- 
leases much of the stress and further expansion of the 
area remarkably leads to further decrease in the stress. 
This "wrong" sign for the extensibility explains why a 
stable pore cannot be achieved under constant tension 
conditions. 

The model of Farago [19] and Tolpekina et al. gives 
perfect agreement with the measured data [57] (see lines 
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in Fig. 12). This, however, not only shows that their the- 
oretical assumptions were correct but also provides a ref- 
erence for comparison of our model with the well known 
model of Goetz and Lipowsky [44, 45] since Tolpekina 
et al. fitted their results to this earlier model. Contrary 
to ours, the model of Goetz and Lipowsky includes ex- 
plicit solvent. That both simulations can be described 
very well by the same theory implies that their physical 
behavior under extensional stress is indeed very similar. 

Fitting the linear extensional regime of our simula- 
tional data to the stretching model we obtain a modulus 
of M w 6.4e/(7^. Translating this to real values depends 
again on the mapping. Using 1.1 e = k^T = 4.1pNnm 
as well as a length mapping of a w 0.9nm gives M w 
30mN/m. This is at the lower end of what one typi- 
cally expects for phospholipid bilayers [21, 23, 42]. The 
rupture tension is Spore ~ 0.55 e/o"^, which translates to 
about 2.5mN/m, while the line tension is 7 ~ 1.2 e/a, 
giving about 5 pN. Both these values are very reasonable 
[46, 47], but one has to be careful, since they depend in 
a nontrivial way on system size [26] . 



VI. CONCLUSIONS 

We have presented a model for the solvent free simu- 
lation of coarse grained lipid bilayer membranes that is 
free of the major problems encountered by earlier efforts 
towards this goal. It robustly self- assembles to a fluid 
phase, uses simple two body potentials, and is highly 
tuneable. We have emphasized that we regard its func- 
tioning to rely on a general principle: the presence of 
sufficiently broad tail attraction potentials. This gives 
lipids a chance of lateral mobility while maintaining flu- 
idity. Since rearrangements and their concomitant local 
increase in pair distances do not immediately cost most of 
the binding energy, the entropy gain upon fluidization is 
not inhibited energetically. In this context we remind the 
reader that it is a well known fact from colloidal physics 
that if the range of attraction is too short compared to 
the particle's hard core radius, no more fluid phase exists 
[48-50]. This principle should therefore also help if exten- 
sions to our present model, such as more beads per lipid 
or a restriction of the long-ranged attraction to speciflc 
beads only are to be employed. 

We have measured many physical characteristics of our 
membrane model and have illustrated that they are easily 
tuneable in a controlled way. Often their values fall well 
inside or close to the experimentally interesting range 
without any explicit careful tuning. This indicates that 
the model is clearly flexible enough and can serve as a 
good starting point for quantitative matching to speciflc 
systems in the spirit of systematic coarse graining, which 
is an approach that has only recently been employed for 
lipid membrane simulations [30, 51] but is well estab- 
lished in other flelds of soft matter [52, 53] . Since we 
have demonstrated that we can readily achieve meso- 
scopic length scales beyond 100 nm and time scales of 



milliseconds, this opens up a wide range of interesting 
mesoscale problems that can with some additional pa- 
rameter matching be simulated quantitatively. Efforts in 
this direction are currently under way. 
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APPENDIX A: SPECTRAL DAMPING 

In Sec. IV A 3 we determined the bending modulus 
from the membrane fluctuation spectrum, which in turn 
was extracted from the simulation by interpolating the 
lipid positions onto a 16 x 16 grid. However, any method 
which interpolates continuous variables onto grid points 
is prone to discretization artifacts, and it is crucial to 
identify them in advance. 

We illustrate the difficulty in the one-dimensional case 
flrst. Let hq{x) be a single mode on the linear interval 
[0;L], given by 

hq{x) = ei(9^+'^<') , g = ^ , n e Z , (Al) 

where ipq is a g-dependent phase. On the grid this mode 
has the values Hq{k) = hq{xk) with Xk = kL/N and k € 
{0,1, . . . , N — 1}. The inverse Finite Fourier Transform 
Hq{u) of these sample points is given by 

f-i 
f-i 

= J_ g27ri(n-«)fe/JVgivPfe 

N ^ 

= e'^'' . (A2) 

This is expected: If n = w we get back the phase of our 
mode, but if n ^ m we get 0, reflecting the fact that the 
wave numbers u and n do not match. 

The situation changes once we go beyond plain sam- 
pling. In the bilayer situation under study our mode 
hq{x) is being represented by many off-lattice points 
X G [0;L]. Grid interpolation was achieved by assigning 
to any grid point the average of the hq {x) for the .t- values 
closest to that grid point. If there are many such points, 
this basically means that we do not sample the mode 
at the grid point, but rather sample its average around 
that grid point. Let us thus deflne the average-sampled 
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interpolation function Hq{k) via 

= ^Qi{QkL/N+v>g) LiqL/2N _ g-ig^,/2JV^ 

iqL V / 

= sinc(^) H,{k) , (A3) 

where we introduced sinc(a;) = sin(7rx)/7r.x. Evidently 
Hq{k) differs from the ordinary sampled interpolation 
function Hq[k) by an additional wave- vector dependent 
but position independent prefactor. This factor goes to 
zero at n = iV, which implies that such modes are en- 
tirely averaged away, but already at the boundary of the 
Brillouin zone, at n 
value sinc^ 

the mesh, it's amplitude is no longer what it used to be 
off-lattice. 

The inverse Finite Fourier Transform of Hq{k) is also 
multiplied by this prefactor: 



Y , it has the disconcertingly small 
^ ~ 0.64. Once the mode has arrived on 



= ^ E sinc(^)//,(fc)e-^-'"'=/^ 

= sinc(^)Hq{u), (A4) 

and thus the power spectrum is damped according to 

\Hq{u)^ = sinc^(^-^^ \Hq{u)\^ = sinc^(^-^) 5n,u ■ 

(A5) 

At the boundary of the Brillouin zone this factor is 
(2/7r)^ « 0.41 and thus not at all negligible. Notice 
that this does not depend on the mode amplitude. In 
other words, even if the mode is only slightly excited and 
the membrane thus looks benignly flat, the estimate for 
the power spectrum derived from the gridded function is 
systematically too low. 

It is easy to see that in two dimensions this result gen- 
eralizes to 

\H,{u)\' = sinc^(f ) sinc2(57) -^-.^ • (^6) 

The fact that this damping is q-dependent renders it 
potentially harmful. Fig. 5 shows that a typical zero- 
tension fluctuation spectrum at low values decays as q~'^ 
but starts to bend upwards once the wave vector ap- 
proaches microscopic scales. The spectral damping dis- 
cussed above will act to suppress this upturn and delay 
it to even larger q-values. The q^^ regime might there- 
fore appear to extend further than it actually does, and 
one might thus be tempted to fit an expression valid only 
in the bending regime to data which owe their exponent 
to a combination of, say, protrusion modes and spectral 
damping artifacts. 



In our mode-analysis we avoided such artifacts by sim- 
ply dividing out the damping factor. 

APPENDIX B: MEMBRANE PORE IN A 
CONSTANT- AREA-ENSEMBLE 

For the total energy of an elastic bilayer spanned inside 

a frame of area A which displays a harmonic extensibility 
with modulus M and which has a line tension 7, [19, 26] 
write 



2 



2\2 



+ 2n'yR 



(Bl) 



where Aq is the tensionless area of that bilayer and R the 
radius of a circular pore. If R is indeed nonzero, its value 
has to be chosen such as to minimize E. 

It is useful to rescale variables. Let us define a char- 
acteristic length scale A and with its help introduce a 
dimensionless pore radius and area extension: 



7^0 
ttM 



R 



R 
A 



and B 



A-Aq 
7rA2 



(B2) 



The equation for the optimal pore radius, dE/dR = 0, 
then reduces to 



R^ - BR+1 = Q 



(B3) 



Hence, the pore opening scenario will exclusively depend 

on only one characteristic dimensionless variable, B; fur- 
thermore, all length scales will be proportional to A with 
a prefactor that depends on B alone. 

It is easy to sec that i? = is always a local minimum 
for E (even though the derivative does not vanish there). 
For B > Be = 3/2^/^ w 1.89 two more stationary radii 
appear as solutions of Eqn. (B3): 



i?±(S) =2/8/3 cos X • (B4) 



The solution R-{B) corresponds to a local minimum 

and for large B asymptotically scales like ^/B; for B > 



= 3/2 



1/3 



2.38 this minimum becomes the global 



one and a discontinuous transition to a pore-state occurs, 
which displays a system size dependent energy barrier of 

£^barrier « 1.38 7A « OM^/'^/^M-^/^AI^^. In this pore 
state the tension is given by 



S_(B) = 1[B- RUB)] ""i' + OiB-^) . (B5) 

At the transition point the tension drops exactly by a 
factor 3, and the system size dependent rupture ten- 
sion is given by Spore = 3/2^3 ^/A rj 2.38 7/A w 
3.49 72/3Mi/3^j^i/3_ 
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